Predicting polarization enhancement in multicomponent ferroelectric superlattices 
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Ab initio calculations are utilized as an input to develop a simple model of polarization in epitaxial short- 
period CaTiOs/SrTi03/BaTi03 superlattices grown on a SrTi03 substrate. The model is then combined with 
a genetic algorithm technique to optimize the arrangement of individual CaTiC>3, SrTiOs and BaTi03 layers 
in a superlattice, predicting structures with the highest possible polarization and a low in-plane lattice constant 
mismatch with the substrate. This modelling procedure can be applied to a wide range of layered perovskite- 
oxide nanostructures providing guidance for experimental development of nanoelectromechanical devices with 
substantially improved polar properties. 



Modern epitaxial thin film techniques make it possible to 
synthesize artificial multicomponent perovskite oxide super- 
lattices whose polar properties can be precisely tailored for a 
wide variety of applicationsji 2 ^ For example, it was recently 
demonstrated that hundreds of atomically thin individual lay- 
ers of CaTi0 3 (CT), SrTiOg (ST) and BaTi0 3 (BT) could 
be grown on a perovskite ST substrate, yielding superlattices 
with compositionally abrupt interfaces and atomically smooth 
surfaces. 2 - 3 It was also shown that — since relaxed lattice con- 
stants of CT and BT are 0.07 A smaller and 0. 1 1 A larger than 
that of ST (ost =3.905 A), respectively — epitaxial strain 
in the constituent layers of these structures can be substantial. 
Due to the strong coupling between strain and polarization 
in ferroelectric perovskites, this can result in substantial en- 
hancement of the polarization relative to that of the bulk con- 
stituents, as has been observed 3 in accordance with theoretical 
predictions^. 

We recently studied such strain-induced polarization 
enhancement in two- and three-component ferroelectric 
CT/ST/BT superlattices epitaxially matched to a cubic 
ST substrate^ First principles methods, namely density- 
functional theory and the modern theory of polarization, 6 
were used to compute the structure and polarization of a small 
number of short-period structures with the same or similar 
compositions as those grown and characterized by Lee et a/. 3 ' 7 
Unfortunately, the substantial computational costs associated 
with these first-principles techniques, growing rapidly as the 
period of the superlattice increases, make it impossible to per- 
form the calculations necessary to answer a broader and more 
interesting question: how should we arrange individual CT, 
ST and BT layers in a given superlattice to obtain the largest 
possible polarization enhancement? In this paper we address 



this question by using our ab initio results as an input to ere- 
ate a simple model for polarization in CT/BT/ST superlattices, 
and then employing this model in conjunction with a genetic 
optimization algorithm technique to identify the optimal can- 
didate structures. 

In previous workjl a simple continuum model was intro- 
duced based on first-principles calculations of ST/BT super- 
lattices; a similar model was subsequently applied to ST/PT 
superlattices. 8 The main premise was to assume that the con- 
stituent layers were linear dielectrics (in the case of the fer- 
roelectric constituent, possessing also a nonzero spontaneous 
polarization), and to obtain the value of the uniform polar- 
ization in each layer by solving the equations of macroscopic 



electrostatics. With appropriate choices for the two dielectric 
constants, this model could reproduce the approximate con- 
stancy of the local polarization in the superlattice, giving a 
nonzero polarization in the ST layer, and the dependence of 
the polarization on the ratio of the thickness of the ST and 
BT layers. However, the electrostatic continuum character 
of this model could not reproduce the dependence of the po- 
larization on the absolute thickness of the constituent layers, 
clearly present in our first-principles results for CT/ST/BT su- 
perlattices. For example, the polarization of the (ST)i(BT)i 
superlattice is noticeably smaller than that of (ST)2(BT)2. 

Here, we introduce a model that includes this "size effect," 
based on the following expression for the energy of the super- 
lattice as a function of the scalar polarization pi of individual 
unit-cell layers i: 



(D 



Here on and f3i describe the anharmonic potential of a sin- 
gle unit-cell layer, and J^j+i represents the coupling between 
nearest-neighborlayers. These parameters take values that de- 
pend on the identity of the layer; for example, a, takes the 
values etc, as, c>b for a CT, ST, or BT layer respectively. 
Similarly, /3j takes the values (3c, (3s, (3b, and there are six 
interface terms J C c, Jbb, Jss, Jsc, Jsb, and J BC . 

To compute the energy for arbitrary values of the unit-cell 
layer polarizations pi would thus require knowledge of twelve 
parameters. However, the approximate constancy of the polar- 
ization across unit-cell layers observed in the first-principles 
results suggests a simplification in which, for each superlat- 
tice, pi is taken to be uniform and equal to the overall polar- 
ization. Substituting p L = p into Eq. Q, we find that 



where 



E{p) = Ap 2 + Bp 4 



A = N u a v + Y N VV ,J VV , 

v («/') 



(2) 



(3) 



(4) 



and N v and N vv i are the number of layers of type v and the 
number of interfaces of type vv' appearing in the superlattice 
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TABLE I: Columns two and three contain expressions for quadratic (^4) and quatric (B) energy-decomposition coefficients for the two- and 
three-component superlattices of Ref. 5. The ab initio ground-state superlattice energies relative to the nonpolar structures, the values of 
ab initio and fitted superlattice polarizations, as well as the differences between the two, are shown in columns four, five, six and seven 
respectively. 



System 


A 


B 


AE (eV) 


|p mi „|(C/m 2 ) 


IiCj < c/ ™ 2 ) 


|Ap min | (%) 


Strained bulk: 














CT 


Jcc 


0C 


-0.019 


0.434 


0.370 


14.85 


BT 


Jbb 


0B 


-0.044 


0.368 


0.363 


1.30 


Two-component: 














(CT)i(ST)! 


2JCS 


0c + 0s 


-0.005 


0.026 


0.026 




/CT^ IRTI 

O l)i (iS l)i 


2 JSB 


PS + PB 


-U.UZJ 


U.ZJ 1 


\J.Zj 1 




(CT)i(BT)i 


2J C B 


0c+0b 


-0.039 


0.231 


0.231 




(CT) 2 (ST) 2 


Jss + Jcc + Ucs 


10c + 2/3 s 


-0.007 


0.168 


0.168 




(ST) 2 (BT) 2 


Jss + Jbb + 2Jss 


20 s + 2I3 B 


-0.039 


0.245 


0.245 




(CT) 2 (BT) 2 


Jcc + Jbb + 2Jcb 


20 c + 20 B 


-0.081 


0.306 


0.306 




Three-component: 














(CT)i(ST)i(BT)i 


JCS + JSB + JCB 


0C+0S + 0B 


-0.034 


0.200 


0.194 


2.98 


(CT) 2 (ST) 2 (BT) 2 


Jcc + Jss + Jbb + Jcs + Jsb + Jcb 


10c + 10s + 20B 


-0.057 


0.242 


0.249 


2.94 


(CT) 2 (ST) 2 (BT) 4 


Jcc + Jss + 3J B b + Jcs + Jsb + Jcb 


2f3 c + 20S + i0B 


-0.131 


0.298 


0.287 


3.83 


(CT) 3 (ST) 3 (BT) 3 


2Jcc + 2Jss' + 2Jbb + Jcs + Jsb + Jcb 


30 c + 30s + 30 B 


-0.082 


0.260 


0.265 


1.91 



sequence. The fact that N c = N C c + [N C s + N C b)/% 
and similarly for TVs an d Nb, for any periodic sequence of 
layers, implies that the three a v parameters and the six J vv t 
parameters enter Eq. (|3j in a linearly dependent way. We can 
then define 

J j a, + aw (5) 

in order to rewrite Eq. l|3} as 

A = N UV J VV , . (6) 

That is, we have eliminated the a„ parameters; from now on, 
we consider our model to be determined by the nine indepen- 
dent parameters (3 C , Ps, Pb, Jcc, Jss, Jbb, Jcs, Jcb, and 
Jsb- 

We obtain the values of the nine model parameters 
{fiv, Juv 1 } by fitting to the first-principles results for the six 
two-component superlattices we considered. For each partic- 
ular superlattice, the quadratic (A) and quatric (B) energy- 
decomposition coefficients, which are the linear combinations 
of {(3 V , J uv i } (see Table[j]for explicit formulas), can be deter- 
mined from first-principles superlattice polarization p m [ n and 
its ground-state energy E(p m [ n ) relative to the structure con- 
strained to have zero polarization. 9 These quantities are re- 
lated to coefficients A and B as follows: 

E( Pmm ) = AE = Ap^ nin + Bp 4 min , (7a) 
=0 A + 2B P 2 min = 0. (7b) 

The resulting parameters </„„'} are shown in Ta- 
blelTTU^ The ab initio ground-state superlattice energies (rel- 
ative to the corresponding nonpolar structures) and polariza- 
tions are shown in the fourth and fifth columns of Table H] 



The values of fitted superlattice polarizations as well as the 
differences between them and their ab initio derived coun- 
terparts are shown in columns six and seven of the same ta- 
ble. The fitted polarization differences for the two-component 
superlattices are not presented, since they are, by construc- 
tion, equal to the ab initio ones. For the rest of the structures 
the model shows a remarkable agreement with first-principles 
results (|Ap m i n | < 4%). For the three-component superlat- 
tices that possess inequivalent polarizations along [001] and 
[001] due to the breaking of inversion symmetry i 5 ^! one could 
in principle compare with the larger polarization, the smaller 
one, or their average. We find empirically that the fit is best 
when compared with the larger polarization, so we have cho- 
sen to present these values in the table. The model performs 
poorly only for the strained bulk CT, whose first-principles 
value of polarization (computed in Ref. 5) represents an ex- 
treme limiting case 12 and cannot be well reproduced by the 
model, which is fitted to the superlattice calculations. 

The availability of such a convenient expression for com- 
puting polarization with nearly ab initio precision allows us 
to predictively identify the arrangements of CT, ST and BT 
layers in a superlattice that would result in the largest pos- 
sible polarization enhancement. While for short period su- 
perlattices {N < 10), this could be done by straightforward 
enumeration, the number of configurations increases rapidly 
with N, necessitating a more sophisticated optimization pro- 
cedure for longer-period superlattices. Here, we use a genetic 
algorithm, 13 in which a particular CT/ST/BT superlattice of 
a given period N is represented by a "chromosome" con- 
taining a sequence of C, S and B "genes". For example, a 
(CT)2(BT)i(ST)2(BT)i superlattice of period 6 is encoded as 
a CCBSSB chromosome. The genetic algorithm also makes 
it easy for us to impose constraints on the optimization, such 
as limiting the thickness of individual layers or the average in- 
plane lattice constant of the superlattice, as discussed further 
below. 
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TABLE H: Fitting parameters used to predict polarization in CT/ST/BT superlattices. See text and Table[I]for more details. 





















Pc 


/3s 


Pb 


Jcc 


Jss 


Jbb 


Jcs 


JCB 


J.5B 


0.584562 


0.648676 


0.833222 


-0.159671 


0.022111 


-0.219844 


-0.000834 


-0.075701 


-0.079061 



Our specific implementation of the genetic optimization al- 
gorithm is as follows. We create an initial population of M 
chromosomes (M is usually in between 2N and 3N) by ran- 
domly assigning C, S or B values to each gene in each chro- 
mosome. The polarization of each chromosome is computed 
using Eq. (7b), and the chromosome's "fitness" is taken to be 
equal to the polarization. The current generation of chromo- 
somes is then replaced by the offspring-chromosome genera- 
tion, created as follows. First, the three chromosomes with the 
highest fitness (the so-called elite chromosomes) are copied 
into the offspring generation without change to preserve the 
best solutions from the previous generation. Second, the re- 
maining M — 3 members of the next generation are created 
by applying the following three-step procedure, (i) Two "par- 
ent" chromosomes are selected from the current generation 
by the so-called "roulette wheel" selection procedure, 14 which 
chooses a chromosome with a probability proportional to its 
fitness, (ii) With probability 10%, the offspring is taken to 
be identical to the parent with better fitness. The remaining 
90% of the time, a "crossover" procedure is applied. We use a 
single-point crossover operator that randomly selects a single 
crossover point on the chromosome and copies the genes from 
one parent up to that point, and from the other after that point, 
(iii) Finally, the offspring is subjected to a "mutation" opera- 
tor, which changes the current value of each gene into one of 
the two other available variants — i.e., gene S, for example, 
could be changed to either C or B — with a low probability 
(in our case: 1%). This entire selection and breeding process 
is continued for five hundred or more generations, after which 
the best available chromosomes are identified. For each set 
of parameters, i.e., the superlattice period and possible layer- 
sequencing restrictions, we perform five separate optimization 
runs to ensure convergence to a consistent solution. 

For any given N, if we impose no restrictions on the num- 
ber of consecutively repeating layers of the same type, then 
the optimal configuration turns out to be pure CT or BT (the 
fitted polarizations of bulk CT and BT are very close, see Ta- 
ble[Q. The former solution dominates in long period superlat- 
tices, while in shorter period ones (N < 10) the latter solu- 
tion is found more often. This happens because, as shown in 



TABLE III: Comparison between ab initio polarizations p m i n and 
fitted polarizations pJJiin f° r a f ew polar short-period superlattices 
identified by the genetic algorithm optimization procedure. 



System 


(JV.fc) 


|p min | (C/m 2 ) 


IPminl (C/m 2 ) 


|Ap min | (%) 


(CT)i<BT) 3 


(4,3) 


0.315 


0.310 


1.6 


(ST) 2 (BT) 3 


(5,3) 


0.279 


0.275 


1.2 


(CT) 2 (BT) 4 


(6,4) 


0.342 


0.328 


4.1 


(CT) 3 (ST)i(BT) a 


(7,3) 


0.313 


0.305 


2.6 



Table U in thin superlattice layers BT has larger polarization 
than CT, which biases the optimization procedure towards BT. 
However, as it is well known, neither of these configurations 
can be experimentally realized because, when grown beyond 
a critical thickness, CT or BT relaxes to its natural in-plane 
lattice constant and the strain-induced polarization enhance- 
ment is lost. Thus, we constrain our optimization procedure 
so that only superlattices containing up to a given number k 
of consecutive layers of the same type are allowed. 

With this "epitaxial growth" constraint, the optimal super- 
lattices that we find fall into two families depending on the 
relation between N and k. For even N and k > N/2 or for 
odd N and k > (N — l)/2, the best solutions have the form of 
(XT) fc (YT)AT_ fc or CXT) fc (YT)^_ fc _i(ST)i, where (X, Y) is 
(B, C) or (C, B). On the other hand, optimal superlattices for 
smaller k (relative to the same period N) contain a number 
of CT/ST/BT stripes and can be reduced to combinations of 
the best solutions of the same form as above but with smaller 
periods. For example, for (N, k) = (12,4) we find three op- 
timal superlattices with polarizations in the range of 0.32- 
0.33 C/m 2 : (CT) 4 (BT) 4 (CT) 2 (BT) 2 , (CT) 3 (BT) 4 (CT) 2 (BT) 3 
and (CT) 2 (BT) 4 (CT) 2 (BT) 4 . Each of these superlattices splits 
into two shorter ones with smaller N and k. These are 
(8,4) (CT) 4 (BT) 4 and (4,2) (CT) 2 (BT) 2 for the first, (7,4) 
(CT) 3 (BT) 4 and (5,3) (CT) 2 (BT) 3 for the second, and two in- 
stances of (6,4) (CT) 2 (BT) 4 for the third optimal superlattice, 
respectively. In what follows we restrict the discussion to so- 
lutions for large k only, assuming that in the opposite case op- 
timal superlattices for any particular N could be constructed 
by merging together an appropriate number of the best large-fc 
solutions for shorter periods. 

We have carried out first-principles calculations for a few 
short-period optimal superlattices to check that their ab initio 
polarizations agree well with those predicted by the model. 
We use a plane-wave based DFT-LDA method 15 with ultra- 
soft pseudopotentiak 16 for structural relaxation of the super- 
lattices and the Berry-phase method of the modern polariza- 
tion theory 6 to compute their total polarization. The details 
of the calculations are the same as in Ref. 5. The results are 
presented in Table ITTTI and show that the good agreement be- 
tween ab initio and fitted values of polarization in short-period 
CT/ST/BT superlattices is preserved. 

Another feature of the superlattice relevant to the feasibil- 
ity of its experimental realization is the mismatch between 
the equilibrium in-plane lattice constant of the superlattice 
(estimated by averaging over the unstrained lattice constants 
of individual layers) and the lattice constant of the ST sub- 
strate. The low substrate mismatch restriction tends to bal- 
ance the number of CT (oct < Qst) and BT (g,bt > 
ost) layers in the superlattice. With this additional screen- 
ing step, we find that the most polar CT/ST/BT superlat- 
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TABLE IV: Short-period superlattices identified by the genetic algo- 
rithm optimization procedure as being the most polar and simulta- 
neously having the lowest lattice constant mismatch |Aa| with the 
substrate. 



System 


i AT l~\ 


iPminl ( L ' m ) 


1 A y-l 1 (G7„\ 

ZACt \70) 


(CT) 3 (BT) 3 


(6,3) 


0.327 


0.34 


(CT) 3 (ST)i(BT) 2 


(6,3) 


0.292 


0.03 


(CT) 3 (ST) 1 (BT) 3 


(7,3) 


0.305 


0.29 


(CT) 4 (BT) 3 


(7,4) 


0.333 


0.12 


(CT) 4 (BT) 4 


(8,4) 


0.337 


0.34 


(CT) 4 (ST)i(BT) 4 


(9,4) 


0.320 


0.30 


(CT) 5 (ST)i(BT) 4 


(10,5) 


0.324 


0.15 


(CT)i (ST)i(BT) 9 


(20,10) 


0.346 


0.24 



tices that emerge from the genetic optimization procedure 
have the following form: (CT) N / 2 (BT) N / 2 f° r even N and 
(CT)( A t_i)/ 2 (ST)i(BT)(jv-i)/2 for °dd N. It is worth point- 
ing out that adding one or two ST layers to CT and BT con- 
taining superlattices destroys their inversion symmetry with- 
out seriously reducing polarization. The lack of the center of 
inversion makes the superlattice polarizations along [001] and 
[001] unequal, which provides for even greater flexibility in 
fine-tuning of the polar properties of such structures. 

In Table IIVI we assemble a number of short-period 
CT/ST/BT superlattices that were identified by the genetic al- 
gorithm optimization procedure as being the most polar su- 
perlattices with a lattice mismatch of less than 0.5%, which 
should allow them to be grown coherently. 7 The following 
first-principles lattice constants were used for the substrate- 
mismatch analysis: acr = 3.813 A (cubic), asr = 3.858 A 
(cubic) and obt = 3.929 A (tetragonal). On average, the po- 
larizations of the superlattices presented in Table HVl are pre- 



dicted to be 10-30% higher than the computed polarizations 
of the previously investigated structures 5 shown in TableU 

To conclude, we have used a first-principles-based one- 
dimensional chain model for polarization in multicomponent 
perovskite-oxide ferroelectric superlattices combined with a 
genetic algorithm optimization procedure to study the con- 
nection between the polar properties of a superlattice and its 
layer sequence. We predict specific layering arrangements 
that produce superlattices simultaneously possessing the high- 
est possible polarization and a low in-plane lattice-constant 
mismatch with the substrate. Our method could be applied 
to superlattices containing individual components other than 
CT, ST and BT, or more than three components, as long as 
the polarization profile across the superlattice remains suffi- 
ciently flat. Various additional restrictions on the arrangement 
of components could easily be added to the genetic algorithm 
optimization to design structures that are custom-tailored for 
specific applications. Our predictions are for ideal structures 
that are defect-free and fully switchable. Since the remanent 
polarization of experimentally grown perovskite-oxide ferro- 
electric superlattices is substantially reduced due to structural 
defects and incomplete switching of ferroelectric domains, the 
computed values are expected to be higher than those ob- 
served. Nevertheless, since our technique does identify the 
most polar layer sequences (regardless of the absolute polar- 
ization) as well as quickly eliminating unfavorable arrange- 
ments, it can be used as a valuable tool to guide the experi- 
mental efforts in the quest for more efficient nanoelectrome- 
chanical devices with tailored and/or substantially enhanced 
properties. 

The authors thank H. N. Lee for sharing his data and many 
valuable discussions. This work was supported by the Cen- 
ter for Piezoelectrics by Design (CPD) under ONR Grant 
N00014-01-1-0365. 
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